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C<l . ABSTRACT 



One of the main results of the Fermi Gamma-Ray Space Telescope is the discovery of 

nj ! 7-ray selected pulsars. The high magnetic field pulsar, PSR J0007+7303 in CTAl, was the 

}nr] I first ever to be discovered through its 7-ray pulsations. Based on analysis of 2 years of LAT 

r-| I survey data, we report on the discovery of 7-ray emission in the off-pulse phase interval at 

Y^', the ~ 6(j level. The flux from this emission in the energy range E > 100 MeV is Fioo = 

P I (1.73±0.40) X 10~* photons cm~^ s~^ and is best fitted by a power law with a photon index of 

(^ I r = 2.54±0.14. The pulsed 7-ray flux in the same energy range is Fioo = (3.95±0.07) x 10"'^ 

photons cm~^ s~^ and is best fitted by an exponentially-cutoff power-law spectrum with a 

photon index of F = 1.41 ± 0.23 and a cutoff energy E^ = 4.04 ± 0.20 GeV. We find no flux 

^ I variability neither at the 2009 May glitch nor in the long term behavior. We model the 7-ray 

i/^ \ light curve with two high-altitude emission models, the outer gap and slot gap, and find that 

the model that best fits the data depends strongly on the assumed origin of the off-pulse 

emission. Both models favor a large angle between the magnetic axis and observer line of 

O I sight, consistent with the nondetection of radio emission being a geometrical effect. Finally 

we discuss how the LAT results bear on the understanding of the cooling of this neutron 

star. 
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Introduction 



The §upernova remnant (SNR) (G119.5+10.2) is a composite SNR characteri z ed by 
a large radio shell enclosing a smaller pulsar wind nebula (PWN). iPineault et al.l fll993[ ) 
derived a kinematic distance of 1.4 ± 0.3 kpc based on associating an H I shell found to the 
northwestern part of the remnant to the remnant itself. Observations in X-rays with ASCA 
and ROSAT revealed a central filled SNR with emission extending to the rad io shell in the 



south and southeast as well as in the north and northwest radio-quiet regions (jSeward et al. 



9951) . X-rav obser vations with ROSAT also revealed the point source RX J0007. 0+7302 



(jSeward et al.lll995l). X-rav observ ations with Chandra revealed a compact PWN and a jet- 



like structure f Halpern et al. 12004). Radio and X-ray characteri s tics o f i mply an age in the 
range 5000-15000 years Jpineault et al.lll993l : Isiane et al.lll997l . l2004h . Isiane et al.l J2004h 
estimated for the age of the SNR a value of 1.3 x 10^(ii,4 yr (where diA is the distance in 
units of 1.4 kpc) which is in g ood agreement with the spin-down age estimate of 14,000 yr 
from 7-rays (jAbdo et al.ll2008l ). The offset of the point source from the geometrical center 
of the radio SNR allows for the estimate of the transverse velocity of the point source which 
Slane et al.l (120041 ) estimated to be ~450 km s~^. 



Prior to the launch of Fermi the EGRET 7-ray source 3EG JOOlO+73 09, which lies 



withi r i the boundaries of t he radio SNR, showed th e characteristics of a pulsar (JBrazier et al. 



19981 : iHalpern et al.l |2004| ) . iMattox et al.l (Il996bl ) discussed this source as a potential can- 
didate for a radio-quiet 7-ray pulsar. A se arch for 7-ray pulsations using EGRET data 
didn't reveal the pulsar (IZiegler et al.ll2008l). The characteris tics of this sou r ce in X-rays 
also pointed to it as a pulsar (IHalpern et al.l |2004| ) . In fact iHalpern et al.l (120041 ) called 
the source the "pulsar" although no pulsati ons were detected from this source using FFT 
searches on ~ 2 6 ks of XMM-Newto n data (jSlane et al.l |2004| ) or in radio with the Green 
Bank Telescope (IHalpern et al.ll2004l ). Very deep searches f or a counterpart for in optical 
and radio resulted only in upper limits (IHalpern et al.ll2004j ). B y correlating X-ray im ages 
from Chandra with those in the optical waveband of the region, IHalpern et al.l (120041 ) gave 
the most accurate position of this source to date ((J2000.0) 00^07""1!56,+73°03'08.1") with 
an accuracy of ~ 0".l. The discovery of the pulsar didn't come until the launch of Fermi. 
During its comm issioning stage th e LAT discovered a 315.87 ms pulsar at the location of 
RX J0007+7303 (JAbdo et al.ll2008l ). The discovery of the pulsar in 7-rays prompted a long, 
~ 130 ks, XMM-Newton observation to search for pulsations from this source (PI: Caraveo 
2008, ObsID: 6049401). Using t his XMM-Newton data set along with the timing model 
from the LAT (Abdo et al.l 2009b) X-ray puls ations from this source in X-rays were finally 



detected ( Lin et al. 



2010 



Caraveo et al.ll2010l ). 



This paper reports further analysis on the LAT data and related observations of the §NR 
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and its 7-ray pulsar, extending the initial report of the pulsar discovery fjAbdo et al.ll2008l ). 
The new developments contribute to a characterization both of the pulsar and its relation 
to the associated extended source, exploiting two years of data now accumulated by Fermi. 
New developments fall in three distinct areas. First, timing analysis of the pulsar spanning 
the two years shows a glitch near the middle of the interval at MJD 54952.652 (1 May 2009), 
with relative change in pulse frequency Au/u ~ 6 x 10~^. This permits careful comparison 
of the pulsar characteristics before and after the glitch. Second, we present new 7-ray results 
concerning the relationship of the pulsar to the surrounding remnant, the first detection of 
an extended source in the off-pulse emission. Finally, we summarize the multi-wavelength 
picture of the source. From this perspective the outstanding characteristic of the neutron 
star in CTAl is that it appears cool for its inferred age. We reconsider both the temperature 
and the age estimates, and t hen relate this to pu lsar characteristics established from timing 
LAT 7-rays. In earlier work ( iHalpern et al.ll2004l ) the pulsar in has been assigned an inferred 
mass exceeding 1.42 Mq even though it is not in a binary, because the larger mass should 
accelerate cooling. Even from initial timing solutions it was clear that this pulsar has a very 
strong magnetic field, though weaker than that of a magnetar, hence it is atypical in two 
respects. 



2. Gamma-Ray Observations and Data Analysis 



We used 2 years of survey data collected with the LAT to study this source in 7-rays. 
The data set starts 2008 August 4 and ends 2010 August 4 (54682.68 - 55412.65 MJD). This 
large data sample w ith high statistics , compared to the six weeks of observations used for 
the dis covery paper (I Ab do et al.ll2008l) and the six months used for the Fermi pulsar catalog 
paper (jAbdo et al.l l2010bl ) allows us to perform several key timing and spectral analyses 
not feasible in prior studies. In particular we study the phase-resolved spectra and the flux 
variability especially around the 2009 May 1 glitch, we search for off-pulse emission and 
build a precise timing model. Throughout this paper we us ed "diffuse" class ph otons events 
with the P6_V11 instrument response functions (IRFs)lll (JAtwood et al.l 120091 ). To reject 
atmospheric 7-rays from the Earth's limb, we selected events with zenith angle < 100°. 



^http://fcrini. gsfc.nasa.gov/ssc/data/analysis/docunientatioii/Cicerone/Cicerone_LAT_IRFs/IRF .overview. html 
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2.1. Timing Analysis 

For the pulse timing analysis, we selected events with energies > 170 MeV that were re- 
const ructed within 1.3° from the Chandra location for RX J0007. 0+7303 source (JHalpern et al 
20041 ) ■ These radius and energy cu ts were selected t o maximize the pulsed significance. Fol- 
lowing the procedure described by iRay et al.l ( 1201 ll ). we measured a total of 72 pulse times 
of arrival (TOAs), each with a integration time of about 10 days, referenced to the geo- 
center. Each TO A was determined using the unbinned maximum likelihood technique from 
Ray et al.l ( 1201 ll ) using a template profile consist ing of 2 gaussian c omponents. The TOAs 
were then fit to a timing model using Tempo2 ( iHobbs et al.l 120061 ). The model, shown in 
Table dl includes position, frequency (z/) and the first three frequency derivatives (z>, i), v'). 
We fitted for position since this will give the smallest timing residuals and will allow for the 
comparison of our timing position to that from Chandra. We find our timing position, shown 
in Table [1], to be consistent with the Chand r a posi tion of compact source RX J0007. 0+7303 
(00:07:01.56, 73:03:08.3; see iHalpern et all J2004h l In addition, a glitch on 2009 May 1 
(MJD 54952.652) was included with an instantaneous step in v and v at the glitch. The 
glitch epoch was chosen to produce a zero phase jump at the glitch . For this glitch we 
measure !S.vjv = 5.54(1) x 10"'^. With a shorter dataset, iRay et al.l (120111 ) could not be 
certain of the step in z> at the glitch. Our extended observations allow us to be confident of 
the Az>/z> = 9.7(6) x 10"^ at the glitch. 

The timing model includes second and third frequency derivatives, which are required 
to obtain white residuals. This is presumably an indication of timing noise in this young 



pulsar. We note that for a braking index of 3, a i> of 1.24 x 10~ s~ is expected from the 
secular spin down of the pulsar. This accounts for about a third of the total measured i). If 
one interprets the measured i) as being entirely due to the secular spin down of the pulsar, 
one obtains a braking index (n) of n = vvji>^ = 9.95. 



2.2. Detection of Off-Pulse Emission 



To search for any 7-ray emission present in the off-pulse part of the phase, we had 
first to determine accurately the definition of the off-pulse phase window. One might simply 
determine the off-pulse interval by eye. In that method the off-pulse starts when any apparent 
pulsed emission decreases to the levels of the background and ends when the pulsed emission 
resumes and an increase above background is seen. A less arbitrary method is to perform a 
likelihood analysis in small phase bins, gradually increasing the width of consecutive bins. 
For all of these width- varying bins only the left and right limits are changing while the center 
of the bin is fixed at what one believes to be the center of the off-pulse interval. A spectral 
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shape is assumed for any off-pulse emission and the increase in signal as a function of phase 
bin width (A0) is analysed. 

We performed this analysis with gtlike where we selected a region of interest of 20° and a 
source region of 30° (see §2.31 for more details). In this analysis all phase bins were centered 
at = 0.89. In the case of absen ce of signal in off-pu lse, one would expect the distribution 
of the test statistic (TS) values (JMattox et al.lll996al ) to be centered around zero with no 
correlation with the width of the phase bins A0. However, in the case of presence of off-pulse 
emission one would expect the significance to increase linearly with the increase in the width 
of the phase bins until one starts integrating photons from the pulsar itself where a sharp 
increase in significance is then seen. The off-pulse width is then defined as the point at which 
the sharp increase in the TS occurs. This is shown in Figure [H 



From the plot, one sees two ranges of data points, with a clear break in slope between 
the earlier range (first seven points) and the remainder. We fitted first order polynomials to 
each range and detected a clear significant change in slope, by about a factor of four. This 
is taken to be the onset of contamination from the pulsed signal. We therefore define the 
off-pulse interval to be G [0.71 — 1.07] with a width of A0 = 0.36. 



As can be seen from Figure [T] there is a significant detection of 7-ray emission in the 
off-pulse phase. At a TS of ~ 40 this is the first detection of 7-ray emission in the off-p ulse of 
PSR J0007-I-7303. An earlier Fermi-LAT survey of PWNe by lAckermann et al.l (1201 ll ). using 
16 months of LAT data, noted a candidate for an off-pulse emission from PSR J0007-I-7303 
but it was below the detection threshold even though a wider phase window was used for the 
off-pulse. The present detection, using 2 years of LAT data but with a conservative phase 
window is unambiguous. 



2.2.1. Off-Pulse Extension Analysis 

Figure [2] shows a TS map of the off-pul se part of PSR J000 7+7303. On the same Figure 
we show ROSAT X-r ay contours in black (jSeward et al.lll995l ) and radio contours in green 
(jPineault et al.l 119971 ). From the Figure one can see that 1) the 7-ray signal detected with 
the LAT is better correlated with the ROSAT X-ray emission than with the radio SNR and 
2) that the off-pulse 7-ray emission is clearly extended. 

To check for a possible extension in the off-pulse 7-ray emission, we performed the likeli- 
hood analysis similar to that in §2.3.1l with the further addition of an extended disk template 
to describe the possible extension of the emission. Different angular sizes in the range 0.1°- 
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1.0° and centroid position templates have been fitted. In each case the significance was 
compared to a simple point-like hypothesis. A disk of radius 0.6° ± 0.3° is favored, and a 
point-like hypothesis is excluded at the 95% confidence level. A template in the shape of an 
eUipse was also fitted to the off-pulse emission. We found no compeUing statistical evidence 
in favor of this fit compared to the disk. 



2.3. Spectral Analysis 

Spectral analyses for this source were performed using the Fermi LAT maximum- 
likelihood Science Tool gtlike in its binned modqj. Fits were performed on a 14° x 14° 
region of the sky centered at the pulsar position selecting photons in the energy range 0.1 
to 300 GeV. We used a model that included diffuse emission co mponents as well a s nearby 



7-ray sources from the First Fermz-LAT 7-ray catalog (IFGL) ( lAbdo et al.ll2010al ) that fell 
within 19° from the position of PSR J0007+7303 . The Galactic diffuse emission was mod- 
eled using the gll_iem_v02_P6_Vll_DIFFUSE model and the isotropic background using the 
isotropic_iem_v02_P6_Vll_DIFFUSE modeE 

In performing the fit we fixed all the parameters of the sources that fell between 14° and 
19° from PSR J0007+7303 to their values in the IFGL catalog, and left free the normalization 
factor of all the sources within 14° of PSR J0007+7303. All the non-pulsar sources were 
modeled with a power law as reported in the IFGL catalog, while the two pulsars in the 
region of interest, PSR J0205+6449 and PSR J2229+6114, were modeled by a power law 
with exponential c utoff according to the data reported in the Fermz-LAT pulsar catalog 



flAbdo et al.ll2010ch . 



To obtain Fermi-LP{T spectral points we divided our sample into logarithmically-spaced 
energy bins (4 bins per decade starting from 100 MeV) and then applied the maximum likeli- 
hood method in each bin. For each energy bin, all point sources, including PSR J0007+7303, 
were modeled by a power law with fixed photon index. From the fit results we then evaluated 
the integral flux in each energy bin. If in an energy bin the source significance is lower than 
3cr we have evaluated the 95% integral flux upper limit in that bin. This method does not 
take into account energy dispersion or correlations among the energy bins. To obtain the 
points of the spectral energy distributions (SEDs) we multiplied the flux in each bin by the 
spectrally weighted mean bin energy. 



^http://fcrmi. gsfc.nasa.gov/ssc/data/analysis/ 

■^ http://fcrnii.gsfc. nasa.gov/ssc/data/acccss/lat/BackgroundModels. html 
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2.3.1. Off-pulse Spectrum 

To quantify the off-pulse 7-ray emission we have assumed a point-hke source, modeled 
with a power law at the pulsar position. We considered only events in the off-pulse phase 
interval [0.71 -1.07]. The fitted power law spectrum is given by: 

dN{E) N{l--i)E-^ 



dE ElnaZ-E^-^ 

where for this fit iV = 1.69±0.40,tat ±0.18,^, x 10"^ photons cm-^ g-^ 7 = 2.54 ±0. 14,^4 ± 
O.OSsys with Emin = 100 MeV and Emax = 100 GeV. The estimated integral flux above 100 
MeV is Fioo = 1.73 ± O.^Ogtat ± O.lSgj^g x 10"^ photons cm~^ s~^ and the integral energy 
flux above 100 MeV is Gioo = 7.83 ± 1.43,^4 ± 0.56,^^, x lO'^^ g^g ^^^-2 g-i_ r^^^^^ -^ ^^ 
compelling statistical case in favor of a cutoff in the spectrum, suggesting that the emission 
is not magnetospheric in origin. We choose to model the off-pulse emission with a power law 
because it has less parameters and gives smaller statistical errors compared to an exponential 
law with a cutoff. 

Systematics are mainly based on uncertainties in the LAT effective area derived from 
the on-o rbit estimations, an d are of < 5% near 1 GeV, 10% below 0.1 GeV and 20% above 



10 GeV (JAbdo et al.ll2009a] ). We therefore propagate these uncertainties using modifled ef- 



fective areas bracketing the nominal ones (P6_V11_DIFFUSE). 



2.3.2. On-pulse Spectrum 

To account for the off-pulse emission we used the results of the off-pulse flt, properly 
rescaled to the on-pulse phase interval, as a starting point for the pulsed emission analysis. 
In the model we considered two sources in the same position, one described as a power law 
with the spectral parameters flxed at the values found with the off-pulse flt and one described 
by a power law with exponential cutoff (PLEC) in the form: 



We set 6 = 1 for which the best-fit parameters are A^o = (9.08 ± 0.20^^* ± 0.54sj,s)x 10 ^^ 
cm-2 s-i MeV-\ T = (1.41 ± 0.23,^4 ± 0.03,^,) and E^ = (4.04 ± 0.20,^^ ± 0.67,^,) GeV. 
The integral flux above 100 MeV is Fioo = (3.95 ± 0.07,^4 ± 0.30,j^,)x 10"^ photons cm-^ 
s~^ and the corresponding energy flux is Gioo = (4.41 ± O.OGstat ± 0.5sys)x 10~^° erg cm~^ 
s~^. We also fitted the on-pulse phase-averaged spectrum to different spectral models, power 
law and broken power law. Both models can be excluded at the > 5a level compared to 
the PLEC model used above. We have also fitted the PLEC spectrum while leaving free 
the exponential index b. The b value obtained was lower then 1 but the overall fit did not 
improve, thus we adopt the simpler PLEC model for which b=l. Figure [3] shows the results 
of the on-pulse phase-averaged spectrum. 



2.3.3. On-pulse Phase-resolved Spectrum 

To explore the on-pulse phase-resolved spectrum, we divided the pulse profile in variable- 
width phase bins, each containing 500 photons above 100 MeV. These bins were defined by se- 
lecting only those events within an energy- dependent radius of 6* < Max(Min(Rmax, 6*68), 0.35° 
around PSR J0007+7303: the minimum value of 0.35° was selected in order to keep all 
high-energy photons, while a maximum radius, Rmax = 2°, was introduced to reduce the 
background contamination at low energies. This choice of binning provides a reasonable 
compromise between the number of photons needed to perform a spectral fit and the length 
of phase intervals. It should be short enough to sample fine details on the light curve, while 
remaining comfortably larger than the rms of the timing solution. A binned maximum like- 
lihood spectral analysis, similar to the analysis performed in §2.3.21 was performed in each 
phase bin with the exception of fixing the spectral parameters of all the nearby 7-ray sources 
and of the two diffuse backgrounds to the values obtained in the phase averaged analysis, 
rescaled for the phase bin width. Using the likelihood ratio test we can reject the power law 
model at a significance level greater than ba in each phase interval. Such a model yields a 
robust fit with a logarithm of the likelihood ratio greater than 150 in each phase interval. 
Figure m shows the evolution of the spectral parameters across PSR J0007+7303's rotational 
phase. In particular, the energy cutoff trend provides a good estimate of the high energy 
emission variation as a function of the pulsar phase. Table [2] summarizes the results of the 
spectral fit in each phase bin. Variations of both the photon index and the cutoff energy 
as a function of rotational phase are apparent. The photon index seems to show a rough 
symmetry centered at a phase point half way between the two peaks with the index increas- 
i ng (softening) out wards. This is similar to what has been observed for PSR J1709— 4429 



(jAbdo et al.ll2010el ). The cutoff energy evolves quite differently as a function of the rotational 



phase. It increases from a minimum value of 1.5 GeV at (/> = 0.2 until reaching the first peak 
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where it stays at a constant value of ~ 4 GeV until the second peak is reached, after which 
it seems to fluctuate between 3 and 6 GeV before finally decreasing to its minimum value of 
1.5 GeV 8X6 = 0.6. 



2.4. Light Curve 

We investigated the pulsar light curve in different energy bands by selecting events within 
1.6° of the pulsar position. This value was selected to maximize the signal-to-background 
ratio over the full energy range {E > 100 MeV). The energy- resolved light curve is shown 
in Figure |3 The top panel in the Figure shows the folded light curve for energies above 
100 MeV while the rest of the panels show the light curve in exclusive energy bands. The 
dashed horizontal line shown in the top panel shows the estimated level of the background 
due to diffuse emission. This background estimate of 194 ± 5 counts/bin was obtained by 
simulating two years of data. We used the LAT Science Tool gtobssim and used for the 
input model the best fitted model from §2.3.21 but with PSR J0007-f7303 and the off-pulse 
component ( §2.3. ip removed from the input model. 

The light curve shows two distinct peaks. We fitted the light curve by a double Gaussian for 
which the first peak and second peak are located at </> = 0.303 ± 0.002 and </> = 0.484 ± 0.002 
respectively. The separation between the means of the two peaks is 0.181 ib 0.003 in phase. 
As can be seen from the Figure there is a significant evolution in the counts ratio of the two 
peaks P1/P2. 



2.5. Flux Variability Analysis 

To check for long-term stability in the flux we performed likelihood analysis similar to 
that in §2.31 but in 8-day time bins. Figure M shows the resulting fluxes. The length of the 
time bin was selected to allow for the accumulation of enough statistics to guarantee a good 
l ikelihood fit. To lo ok for flux variability from the source we adopt the method outlined in 



( lAbdo et al.ll2010al ). The source shows no variability on this time scale. Similar analyses 



were performed for 16, 32, and 64-day time bins, no significant modulation was found. 

To check for any change in the spectrum of the pulsar due to the glitch we split the 
data in two bins around the glitch. The pre-glitch epoch spans the time range 54682.68 
- 54952.652 (MJD), while the post-glitch data spans the time range 54952.652 - 55412.65 
(MJD). We performed a likelihood analysis similar to that in §2.5l in these two time bins. In 
Table E] we show the spectral fits for the pulsar in these two epochs. The flux and spectral 



-10- 



parameters are in good agreement for the two epochs. No change in integral flux above 100 
MeV is seen. 

Recent variab ihty detected in the Crab appears to come from the Nebula, not the pulsar 
(JAbdo et al.ll201ll ). Since PSR J0007+7303 is among the younger pulsars and has a complex 
PWN/SNR associated with it, it is reasonable to ask whether analogous variability occurs 
in this source. While §2.21 has shown strong evidence for an off-pulse component, tentatively 
extended, there are not enough statistics to explore variability in this component at all, and 
certainly not on the time scales detected in the Crab. 



3. Discussion 



3.1. Geometrical Constraints from Light Curve Modeling 

The pulsar emission mechanism is not well understood, and the distributions of pulsar 
magnetic inclination angles, efficiencies, and other physical characteristics are largely un- 
known. There are several competing emission models, of which the resulting pulse profiles 
depend on the geometry of the system, defined by the emission zone location and size, ob- 
server viewing angle, and magnetic inclination angle. Fitting the profiles of these models 
to observed light curves can provide insight on the true emission and viewing geometries, 
for example leading to better constraints o n luminosity and effic iency through calculation of 
the fiux correction factor /q (equation 4 of IWatters et al.l (120091 )). The geometry of a given 
system likely also contributes to the detectable presence or absence of radio and perhaps 
X-ray emission. For example, if the inclination and observer angles are very different, we 
could see a 7-ray-only pulsar, as the narrower radio beam may not cross our line of sight. 
To determine the geometry of the CTA 1 pulsar and distinguish between emission models, 
we compared the LAT light curve of PSR J0007-I-7303 with the predicted light curves from 
geometrica l representations of two stand ard high-energy pulsar emission models, the outer 
gap (OG) (IRomani fc Yadigaroglulll995[) and slot gap (SG, also referred to as the two-pole 
caustic or TPC) (JMuslimov fc Harding] |2004| ) models. These models were considered within 
the context of the vacuum retarded dipole magnetic field. 



The light curves were simulated as in lDyks et al.l (120041 ). with the geometry modified to 
represent each emission region. The OG lies along the last open field lines between the null 
charge surface and the light cylinder, R\c = c/Q, with rcyi < 1.0i?ic as large as possible for 
each geometrical configuration. The SG extends from the surface to r^ji = 0.95i?ic. In both 
geometries, the maximum emission altitude r^ax was treated as a model parameter. We note 
that our representation of the SG differs from that of the TPC model in lRomani Sz Watters 
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( 120101 ) , as we allow emission out to much higher altitudes (their emission is cut off at rcyi = 
0.75-Ric). The TPC geometry is therefore treated as a subset of the SG, rather than being 
considered as a completely separate model. The SG model has emission throughout the gap, 
while the OG emission occurs only along the field line at the gap's innermost edge. In the 
simulations, the vacuum retarded dipole field is assumed in the observer's frame and then 
transformed to the co- rotating frame (CF). Photons are emitted tangent to the field in the 
CF, prior to the aberration calculation. A constant emissivity is assumed along the field 
lines in the CF. The special relativistic aberration leads to a bunching in pulse phase, or 
caustics, on the trailing side of the pulse, producing peaks in the li ght curvyes . For a given 



inclination angle a, gap width w in open volume units {rove, as in iDyks et al.l (J2004f )). and 
maximum emission radius rmax in units of Ric, the code produces light curves at all observer 
angles (. 

The LAT light curve of PSR J0007+7303 has 32 bins and a background count level 
of 195 or 246 counts/bin, depending on whether the off-peak emission discussed in §2.2 is 
(respectively) magnetospheric or due to a PWN. We compared this light curve with light 
curves simulated from a representative sample of all possible geometries and rebinned to 
match the observed light curve. Our models had five free parameters, a, (, w, rmax, and a 
phase shift A0 introduced in order to best match the LAT light curve. We considered values 
of a between 0° and 90° and ( between 0.5° and 89.5°, each with resolution of 1°; values of 
width < u; < 0.3 with resolution 0.01 TovcJ and maximum emission radii 0.7 < rmax ^ 2 
for the SG and 0.9 < r < 2 for the OG, with resolution 0.1 -Ric in both cases. For each 
model, the emission altitude is limited by Min(rinax5 '"cyi) — '"max therefore may be larger than 
rcyi, but emission ceases at the cylindrical radius if not at the maximum emission radius. 
The phase shift 0° < A0 < 360° is arbitrary because PSR J0007+7303 is radio-quiet within 
the flux limits achieved thus far by radio telescopes; were the pulsar radio-loud, the shift 
would be constrai ned to be at most equal to the phase lag between the radio and 7-ray peaks 



(IDvks et al.ll2004h . 



To find the best-fit parameters of each model, w e used a Ma r kov C hain Monte Carlo 



(MCMC) maximum likelihood routine as described in I Verde et al.l ( 120031 ) that explored the 
parameter space. We calculated the x^ between the LAT light curve and model light curve 
and used Wilks' Theorem, A ln(L) = — Ax^/2, to guide the Markov chains toward regions of 
high likelihood. We also used the Ax^ test to find 3a confidence intervals on each parameter. 
For this pulsar, the best-fit light curve has a very large x^- This is expected, as a) the pulsar 
is bright and its light curve has relatively small error bars and b) the light curve simulations 
result from geometrical models of possible but simplified pulsar magnetospheres, rather than 
from a well-understood physical model. When finding confidence intervals in steps of Ax^ 
from a large initial x^ value, the intervals will appear artificially small. We therefore re-scale 
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the x^ values found for all sets of parameters in the MCMC chains so that the minimum 
reduced y^ = A^dof = 27 (the scale factor is therefore N^oi/Xmm)- The confidence intervals 
are then found with these modified x^ values. In this way, we find parameter ranges that 
may give rise to the high-energy light curve of PSR J0007+7303. 

For the case where the off-peak emission is assumed to be magnetospheric (corresponding 
to a background level of 195 counts/bin), the best-fit parameters are (a, (, w, r, A0) = (6°, 
74.5°, 0.04rovc, 1.0i?ic, 2°) with xV27 = 22.5 for the OG model and (8°, 69.5°, Otovc, 1.9i?ic, 
-4°) with x^/27 = 11.3 for the SG model. Almost identical parameters are found for the case 
where the off-peak emission is assumed to originate from a PWN. The confidence intervals, Jq 
values, and reduced x^ for the best fit in each interval are given in Table HI The absolute best 
fit is found with the OG model using the higher background level. Figure [7] shows the LAT 
light curve superposed with the best-fit model light curves, as well as reduced x^ contours in a 
and ( with w and r fixed at the best-fit values. To be complete, we also fit the light curve with 
the TPC model, the results of which are not shown here; we find this lower-altitude subset of 
the SG model cannot reproduce the sharp double peaks with the correct spacing anywhere 
in parameter space. We therefore find that both outer magnetosphere models produce light 
curves consistent with that which is observed. For magnetospheric off-peak emission, the SG 
is preferred over the OG; its ability to reproduce the light curve is seen clearly by comparing 
the red curves in panels (a) and (c) of Figure [3 The reverse is true under the assumption 
of off-peak emission from a wind nebula, shown by the blue curves. Under the assumption 
of magnetospheric off-pulse emission, the SG misses some emission in the wings and has too 
high a background level; this effect is magnified when we assume instead a PWN origin of this 
emission. Such details may be modified with a more physical emission model, for example 
by including azimuthal asymmetry in the accelerati ng electric field from offset polar caps. 



which leads to a decrease in the off-peak emission (JHarding &: Muslimovll201ll ). The OG 
characteristically does not reproduce the wings or higher emission in the off-peak phases 
when the background is derived assuming pulsar emission in the off-peak; by definition it 
matches the background well when the background is instead found assuming the emission 
is not from the pulsar itself. 

Regardless of the background counts used, the best-fit values of a, and ( for the OG and 
SG models are far apart, with |C~Q;| >60 degrees. The radio beam would need a width of that 
order in order for the radio pulse to be detect able. For the para meters of PSR J0007-I-7303 a 



model estimate of beam width is < 10 degrees (jStory et al.ll2007l ). Thus the preference for OG 



and SG model fits over that of TPC, along with the fitted parameter values, self-consistently 
offers a satisfactory explan ation for non-detection of any radio pulse. A deep search using 



GBT (JHalpern et al.ll2004l ) yielded an upper limit that remains the lowest limiting flux for 
any pulsar position. The upper limit on the pseudo-luminosity of L1400 — 0.02 mJy kpc^ is 
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lower than the lowest pseu do-luininosity measu red for a radio pulsar (L1400 — 0.035 niJy 
kpc^ for PSR J1907+0602 fJAbdo et al.ll2010dl )) and is thus restrictive enough to support 
the radio-quiet designation. Conversely, if one were to adopt the TPC model for this source 
then the nearly-equal values of a and C would suggest that a radio pulse should have been 
found. 

X-ray pulsations were recently detected from the CTA 1 p ulsar using the same LAT 
ephem eris used in the timing and light curve analysis of this paper ( JLin et al.ll2010t ICaraveo et al, 
2OIOI ). We did not include any X-ray information in our fits, but can consider whether or not 
our results make sense in the multiwavelength picture. The peak in the X-ray light curve 
shows a significant thermal component, suggestive of a hot spot on the surface. The authors 
estimate the size of the hot spot to be ~ 100 m in radius, larger than the polar cap in the 
case of a dipole model for the pulsar; it is not clear from the data where the hot spot is 
located. The X-ray peak occurs at ~ 0.25-0.3 (~ 90°-110°) prior to the first peak in the 
7-ray light curve. Our HE models predict a similar location of the magnetic pole: the first 
7-ray peak is at 90°, and in all our model fits the pole lies 36°-94° before this peak; taking 
only the SG and best OG fits places the pole 94°-88°. Our best-fit models are therefore 
consistent with the thermal X-ray emission originating near the magnetic pole. 



3.2. Thermal Component in X-rays and Neutron Star Cooling 

PSR J0007-I-7303 is among the neutron stars (NS) where the theory of cooling confronts 
observation. It held that distinction before launch of Fermi, based on the candidate NS 



identified at that time (IHalpern et al.ll20p4l). a nd it continues to be treated as an exceptional 
case in the cooling literature (jPage et al.ll2009l ). Those references cite additional literature on 
heat loss models which describe how observable surface temperature or thermal flux declines 
with age. X-ray or extreme ultra violet surface emission in young, hot NS provides the 
observational test. The NS age and surface temperature must be known or constrained, that 
is, a limit may be useful if corrections can be characterized at least as to sign if not magnitude. 
Only comparatively young NS within a few kpc provide useful constraints. The only pulsar 
both younger and closer than PSR J0007+7303 is the Vela pulsar. Before Fermi, the age 
of PSR J0007+7303 was equated to the estimated CTAl SNR age and th ermal flux was 
estimated from the X-ray emission of the identified point-like NS candidate (IHalpern et al. 
20041 ). From this it was understood that the putative pulsar in CTAl was cool for its age, in 
comparison with both theory and other young NS. This led to theoretical effort to understand 
why this particular NS had cooled rapidly. 



Fermi results impact context information that enters into the analysis of the cooling 
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question. The pulsar ephemeris affects the X-ray analysis used to extract the thermal ffux. 
The Fermi results have not greatly altered best-estimates of key quantities but provide a new 
network of interlocking constraints that give greater robustness of the observational context 
information, as follows. (1) The spin period P and its derivative P give spindown age of 
~13,900 yr for the pulsar. This age estimate is well within the broader range of historical 
estimates for the SN R age, 5,000-15,0 00 yr, and nearly identical to one pre-launch estimate 



of 1.3 X 10^(ii4 yr (jSlane et al.ll2004j ). The formal propagated error in the spindown age 
estimate is negligible. The principal issue is validity of the dipole approximation formula 
but it clearly reinforces the age estimates derived other ways. (2) The same observations (P, 
P) also yield spindown energy loss independently of distance. There are no inconsistencies 
between this energy budget and that of the pulsar plus SNR, using the accepted distance, 
D, of 1.4 ± 0.3 kpc, hence the pie-Fermi distance estimate used is not in conflict with new 
information and the distance scale factor in the age estimate just cited is consistent with 
unity. (3) The pulsar ephemeris has now been used to detect X-ray pulsations using XMM 



( jLin et al.ll2010l : ICaraveo et al.l 120101 ). The X-ray spectrum can be divided into pulsed and 
unpulsed components. While total X-ray flux is consistent with levels measured in pre- 
Fermi work, the part now assigned to a DC thermal component emitted from the surface is 
reduced by at least a factor that equates to the unpulsed fraction. Since D is unchanged, the 
revised bound on NS surface temperature remains as in earlier estimates or perhaps even 
slightly reduced. One may consider modiflcations due to the interstellar absorbing column 
and to the neutron star atmosphere causing the surface flux to depart from blackbody. 
These considerations properly enter into many other NS used for comparison with cooling 
and it is beyond scope of this paper to reanalyze the entire X-ray discussion. It is left to 
future research but provisionally the pre-launch thermal luminosity stands. (4) In this same 
connection a further point to note is that the off-pulse 7-ray flux, discovered and quantifled 
in §2.21 of the present paper, cannot be thermal flux from the neutron star and must be 
taken as magnetospheric or PWN emission. In either case, that same magnetosphere or 
PWN component could extend downward in energy to X-rays, and in principle could provide 
further downward adjustment to the thermal X-ray flux from the star. This is left as an 
adjustment of unknown magnitude but known sign; it can only require the star to cool faster 
than the estimate obtained by neglecting it. (5) From P and P one also obtains an estimate 
of the stellar dipole fleld, B ~ lO^^G. PSR J0007+7303 therefore falls among the most highly 
magnetic neutron stars, magnetars excepted. Cooling scenarios whereby a strong magnetic 
fleld could mod ify the cooling curve h ave been described in earlier literature (for a summary 



see for example lYakovlev et al.l ( 1200 ll )). The high value of B now established for this pulsar 
means mechanisms whereby high B enhances cooling may merit further attention. (6) The 
previous section shows how model flts for OG and SG models favor small a and large C, and 
could explain the absence of radio pulsations. 
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The first four items in tlie list mean tliat if PSR J0007+7303 was an outlier relative 
to models before Fermi and relative to pTe-Fermi theoretical understanding, it has become 
slightly more egregious relative to that prior theoretical understanding. However theory has 



also advanced, largely from observation of cooling in the central star in Cas A (jPage et al. 



20091 ). Data from that source are now the most definitive and constraining of any un- 
recycled NS. Fitting Cas A has given prominence to NS cooling models involving Cooper 
pairing contributions. Points (5) and (6) in the list regarding magnetic field strength and 
geometr y may provide guida nce to theory. Sufficiently strong magnetic fields can affect 



cooling ( lYakovlev et al.ll200ll ) and might affect heat flow in the star. Also the particular 
geometry in PSR J0007+7303 with low a and high ( that operates against radio pulse 
detection could also lead to a misleadingly low thermal X-ray flux if the strong field conveyed 
internal heat flow preferentially to the magnetic polar regions. Then, a Lambertian emission 
pattern from the hotter poles would be anisotropically beamed away from an observer at 
high (. Further X-ray observations, combined with multi-wavelength analysis applied and 
field geometry modeling may shed further light on the thermal luminosity. This could be 
undertaken comparatively with other young, nearby pulsars such as the Dragonfly pulsar, 
PSR J2021+3651, which has estimated age 17 kyr, distance 2.1 kpc, dipole field 3 x 10^^ 
G, and has radio pulses. Such comparisons might bring out the role of the magnetic field 
strength and geometry in NS cooling. 



4. Summary 



PSR J0007+7303 is among the brightest 7-ray pulsars (jAbdo et al.ll2010d ). It also is 
part of an interesting PWN and SNR complex that is still young. We have exploited greatly- 
improved cumulative statistics from two years of Fermi-LAT data to investigate questions 
that can be pursued only on brighter 7-ray pulsars. Interesting aspects of the source include 
its having had a major glitch during the Fermi observing period, its being among the most 
strongly magnetic of neutron stars that are not magnetars, and its status as a prime testbed 
for neutron star cooling theory, a topic now receiving renewed attention because of the 
cooling observed in X-rays in Cas A. 

Pulsar phase dependence of the light curve and spectrum have been investigated and the 
source has been compared both with other well-studied pulsars (notably PSR J1709-4429) 
and also with standard magnetospheric geometry models, revealing a clear preference for 
models where emission occurs high in the magnetosphere; particularly the SO model. Glitch 
parameters have been extracted. We have conducted a systematic search for long-term vari- 
ability in the system, with negative results. Neither a change associated with the glitch 
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nor flaring such as has recently been seen in the younger Crab Nebula f lAbdo et al.ll201l[ ) 
has been found. However, off-pulse emission has finally been detected at high confidence 
and there is evidence that it is extended; hence there is now a new 7-ray component in 
the overall source, potentially a PWN although the possibility that it originates inside the 
magnetosphere is not strongly excluded. The variability of that off-pulse source is a subject 
to be pursued as Fermi continues to accumulate data. We have described how the parame- 
ters emerging from the 7-ray analysis (ephemeris, spindown energy loss, age, distance, and 
magnetic field) and from follow-on X-ray studies affect understanding of the cooling history, 
reinforcing the conclusion that the surface of this NS is cool for its age. This now needs to be 
followed up with additional X-ray analysis to further constrain the surface X-ray luminosity. 
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Fig. 1. — Test Statistic trend versus phase-bin widths. The dashed hne represent the detec- 
tion threshold (TS = 25). In all of the bins the center value was = 0.89. 
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Fig. 2.— Fermi-LAT TS map of the off-pulse part of PSR J0007+7303 (0 G [0.71 - 1.07]). 
Black cross marks t he location of the pulsar. ROSAT X-ray contours are sh own in black 
(ISeward et al.lll995l ). Radio contours are shown in green (jPineault et al.lll997[ ). 
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Fit and data-set 



Pulsar name J0007+7303 

MJD range 54682.7—55415.4 

Number of TOAs 72 

Rms timing residual (ms) 2.3 

Reduced x^ value 1.39 



Measured Quantities 



Right ascension, a 00:07:01.7(2) 

Declination, 5 +73:03:07.4(8) 

Pulse frequency, v (s"^) 3.165827392(3) 

First derivative of pulse frequency, z> (s~^) — 3. 6120(5) x 10"^^ 

Second derivative of pulse frequency, D (s~^) 4.1(7) xlO~^^ 

■z>' (s-4) 5.4(9) xlO-30 

Az/ 0.000001753(2) 

Az> -3.5(2) X 10-15 



Set Quantities 



Epoch of frequency determination (MJD) 54952 

Epoch of position determination (MJD) 54952 

Epoch of dispersion measure determination (MJD) 54952 

Glitch Epoch 54952.652 

Phase jump at glitch 



Derived Quantities 


logio (Characteristic ^&^i Y^) 

log j^o (Surface magnetic field strength, G) 


. 4.14 
. . 13.03 


Assumptions 


Solar system ephemeris model 


. . DE405 



Table 1: Measured and Derived timing parameters of PSR J0007+7303. Figures in paren- 
theses are twice the nominal la tempo2 uncertainties in the least-significant digits quoted. 
The time system uses is TDB. 
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Fig. 3. — On-pulse phase-averaged spectral energy distribution of PSR J0007-I-7303. The 
sohd black line represents the best fit power law with exponential cutoff with b=l. Dashed 
lines represent the la errors in each case. 
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Fig. 4. — Evolution of photon index (top) and energy cutoff [bottom) above 0.1 GeV as a 
function of pulse phase from fits in fixed-count phase bins of 500 photons per bin. The error 
bars denote statistical errors. For each phase interval (defined in Table |2]) a power law with 
exponential cutoff has been assumed. The dashed histogram represents the LAT light curve 
above 0.1 GeV. 
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Fig. 5. — The evolution of the pulse profile of PSR J0007+7303 with energy. Two rotational 
periods are shown with a resolution of 32 phase bins per period. The top panel shows folded 
light curve for energies above 100 MeV. The dashed horizontal line shown in the top panel 
shows the estimated level of the background due to diffuse emission (see text for details). 
The rest of the panels show the light curve in exclusive energy bands. The darker histogram 
on the second panel from the top shows the folded light curve for energies > 5 GeV 
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Fig. 6.— Flux (> 100 MeV) of PSR J0007+7303 as a function of time in 8-day time bins. 
The flux shows no evidence for variability. The dashed vertical line marks the time of the 
ghtch. 
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Fig. 7. — The best fitting model light curves for each model, red indicating 195 counts/bin 
and green 246 counts/bin as the background level, and the reduced Xr = id I^Aai contours 
in a-C, parameter space for the case of a background of 195 counts/bin. (a) Best fitting 
outer gap light curve (red/blue) superposed on the LAT light curve of the PSR J0007+7303 
pulsar (black); the best fit parameters are given in the text. The horizontal dashed lines 
represent the estimated background levels, while the vertical dotted lines mark the location 
of the magnetic pole in the context of the outer gap model for the listed sets of parameters. 
(6) Xr contours for outer gap light curves in which w and r have been fixed at the best-fit 
parameters. The filled contours show x^ in increments of 20 for all a and C- The best fit is 
marked by the pink triangle. The 3 a confidence intervals for each parameter are given in 
Table m (c) Same as (a), for the slot gap model, (c?) Same as (&), for the slot gap model. 
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(Pmin 


(Pmax 


Photon index 


Cutoff energy 
(GeV) 


Flux (> 100 MeV) 
(x 10"'''photons cm"^ s~^) 


0.07 
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0.882 ± 2.419 
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1.264 ± 0.091 


4.116 ± 0.779 


9.575 ± 0.91 


0.443 


0.457 


1.333 ± 0.085 


4.897 ± 0.966 


10.91 ± 1.003 


0.457 


0.471 


1.249 ± 0.085 


4.755 ± 0.866 


10.12 ± 0.9369 


0.471 


0.486 


1.321 ± 0.078 


5.788 ± 1.124 


10.768 ± 0.9581 


0.486 


0.5 


1.111 ± 0.096 


3.178 ± 0.519 


10.142 ± 0.95 


0.5 


0.513 


1.161 ± 0.096 


3.095 ± 0.523 


11.402 ± 1.042 


0.513 


0.528 


1.398 ± 0.081 


5.504 ± 1.164 


10.126 ± 0.926 


0.528 


0.547 


1.483 ± 0.099 


3.449 ± 0.772 


9.344 ± 0.892 


0.547 


0.574 


1.602 ± 0.121 


3.218 ± 0.921 


5.652 ± 0.649 


0.574 


0.623 


1.48 ± 0.263 


1.502 ± 0.595 


1.635 ± 0.356 


0.623 


0.71 


0.762 ± 2.746 


0.689 ± 1.44 


0.332 ± 0.405 



Table 2: Phase interval definitions and corresponding spectral parameters obtained from 
fitting the spectrum with a power law with exponential cutoff. The flux in the third column 
is normalized to the width of the phase bin. The systematic uncertainties are in agreement 
with the ones evaluated for the phase-averaged analysis. 
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Pre-glitch Post-glitch 



Date range (MJD) 54682.68 - 54952.652 54952.652 - 55412.65 
Photon index 1.42 ± 0.04 ± 0.03 1.50 ± 0.03 ± 0.04 

Cutoff energy (GeV) 4.65 ± 0.39 ± 0.77 4.74 ± 0.28 ± 0.79 

Flux (> 100 Me V) 3.60 ± 0.11 ± 0.27 4.09 ± 0.08 ± 0.31 



Table 3: Phase-averaged spectral parameter and flux for PSR J0007+7303 for the two epochs 
around the glitch. Flux is given in 10~^ photons cm~^ s~^ units. First errors are statistical 
and second ones are systematical errors. 



Parameter 


Outer Gap 1 


Slot Gap 1 


Outer Gap 2 


Slot Gap 2 


«n 


6+2> 83+^ 


8+^ 


Qtl 


84, 8 


en 


74.5+^^ 13.5+^1 


69.5+° 


79.5tl 


11.5l^, 69.5 


w (rove) 


OMt'o-^l omtr 


o.oir' 


09+°-°^ 


0.03+r'> 


r (Ric) 


i.otli, i.otl-' 


> 1.0 


1 1+0.9 
-'-•-'--0.1 


0.7, > 1.0 


A0n 


2tl 5it% 


_4+4 


1Q + 10 


186+g, 356 


XV2T 


22.5, 24.5 


11.3 


7.0 


9.3, 21.9 


h 


0.17, 1.10 


0.42 


0.15 


2.0, 0.42 



Table 4: Best-fit model parameters of the OG and SG geometrical models to the LAT light 
curve of PSR J0007+7303, where the background used for the fitting was 195 counts/bin for 
columns labeled "1" and 246 counts/bin for those labeled "2". The asymmetric error bars 
give the 3a confidence intervals, derived using the scaled x^ as described in §3.11 There were 
two regions in parameter space that fell within 3a of the best OGl and SG2 fits; the first 
set of values given in the column are the best of the two regions. The quoted reduced x^ 
and /q correspond to the absolute best fits in each region. The parameters shown here are 
used in Figure [7] for the simulated light curves and x^ contours. 



